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Abstract 

We present the package SADE (Symmetry Analysis of Differential Equations) for the 
determination of symmetries and related properties of systems of differential equa- 
tions. The main methods implemented are: Lie, non classical, Lie-Backlund and 
potential symmetries, invariant solutions, first-integrals, Nother theorem for both 
discrete and continuous systems, solution of ordinary differential equations, order 
and dimension reductions using Lie symmetries, classification of differential equa- 
tions, Casimir invariants, and the quasi-polynomial formalism for ODE's (previously 
implemented by the authors in the package QPSI) for the determination of quasi- 
polynomial first-integrals, Lie symmetries and invariant surfaces. Examples of use 
of the package are given. 
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PROGRAM SUMMARY 



Title of the program: SADE 

Catalogue identifier: None 

Program obtainable from: The author by e-mail. 

Operating systems under which the program has been tested: UNIX/LINUX 
systems and WINDOWS 

Programming language used: MAPLE 13 and MAPLE 14 

No. of bytes in a word: 32 

No. of bytes in distributed program: 300 KB 

Distribution format: zip or gzip 

Card punching code: ASCII 

Keywords: symmetry transformations, invariant solutions, first integrals, nother 
theorem. 

Nature of the physical problem: Determination of analytical properties of sys- 
tems of differential equations, including symmetry transformations, analytical 
solutions and conservation laws. 

Method of resolution: The package implements in MAPLE some algorithms 
(discussed in the text) for the study of systems of differential equations. 

Restrictions on the complexity of the problem: Depends strongly on the system 
and on the algorithm required. Typical restrictions are related to the solution 
of a large over-determined system of linear or non-linear differential equations. 

Typical running time: Depends strongly on the order, the complexity of the 
differential system and the object computed. Ranges from seconds to hours. 

LONG WRITE-UP 

1 Introduction 

Natural phenomena are very often modeled by differential equations, which 
exhibit a plethora of dynamical behaviors. These can be classified somewhat 
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vaguely in two categories: regular and irregular, according to the complex- 
ity exhibited by its solutions. The notion of integrability is then used as an 
attempt to put a more stringent distinction between regular and irregular 
behavior. Usually a regular behavior is characterized by the existence of con- 
servation laws that strongly restrict the types of solution a system can exhibit. 
Even for non-integrable systems some of these laws can be obtained. Also the 
determination of particular analytical solutions for both Ordinary (ODE's) 
and Partial Differential Equations (PDE's) is of utmost importance in many 
fields of physics and applied sciences. These solutions are helpful to shed some 
light and gain insight on the physics of the system, and are also useful as 
benchmarks for numerical methods. Almost all known analytical solutions in 
physics are solutions invariant under one or more symmetry transformations. 

The theory of symmetry transformations of Differential Equations (DE's) was 
introduced by Lie in the end of the XIX century [I]. Solutions invariant under 
symmetry transformations are called invariant solutions, and different meth- 
ods are described in the literature (see References [2, 3,4,5] and references 
therein). Olver and Rosenau showed that new solutions can be obtained by re- 
quiring that they are invariant under infinitesimal symmetry transformations 
while also preserving additional side conditions [3J. Non classical symmetries 
were introduced by Bluman and Cole and are based on the idea that the 
required analytical solution is invariant under symmetry transformations pre- 
serving both the form of the differential equation and the invariant solution 
condition [5J . This approach is less restrictive in the sense that there exist usu- 
ally more non classical symmetries than Lie (classical) symmetries, the latter 
being a subset of the former. Other generalizations of the classical Lie method 
considered here are potential symmetries [6] and Lie-Backlund transforma- 
tions [7j. For a first introduction to Lie Symmetries see [8f9] . and pT0][2] for 
a more complete and formal approach. A description of methods for solving 
differential equations using Symmetries is found in 

Different packages in computer algebra systems exist implementing Lie sym- 
metry computations and related methods: SPDE by Schwarz [12], CRACK 
and LIEPDE by Wolf [131114] and DIMSYM by Sherring and Prince [15] in 
REDUCE, LIE and BIGLIE by Head [TdUTT] in MUMATH and MATHLIE by 
Baumann [IB] in Mathematica. For MAPLE there are also some useful pack- 
ages: PDEtools by Cheb-Terrab [19] which is distributed since Release 11, 
DESOLV by Vu and Carminati [20114] . and GeM by Cheviakov [21]. For good 
reviews with a comparison between some of these packages see References [22] 
and [23] . The package QPSI by the authors implements the Quasi-Polynomial 
formalism for symmetry generators, first-integrals and invariant hyper-surfaces 
for ODE's [24], now part of the present package. 

In this work we present the package Symmetry Analysis of Differential Equa- 
tions (SADE) in MAPLE, for the computation of Lie, Lie-Backlund and non- 
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classical symmetries, invariant solutions, first-integrals, Nother theorem for 
both discrete and continuous systems, quasi-polynomial first-integrals and 
symmetry generators, solution and reduction of order or dimension for ODE's, 
classification of differential equations, invariant surfaces and Casimir invari- 
ants [8,10,25,24], and some other features presented below. Our package is 




well suited for efficiently computing Lie symmetries of large systems, as for 
instance the Yang-Mills with SU(2) and SU (3) gauge group [26J, and has been 
used in the last years in our group in different applications [27,28, 29, 30l3"T]. 

Our aim in developing this package was to implement these methods, being 
as user friendly as possible for researchers in many fields of pure and applied 
sciences, and still being capable to handle reasonably complicated systems of 
equations. The paper is structure in the following way: in section 2 we briefly 
revise the mathematical methods implemented. A discussion of the heuristics 
for the solution of the determining system for Lie and nonclassical symmetries 
of linear and non-linear overdetermined systems of PDE's is given in section 3. 
The package routines are described in section 4, and some illustrative examples 
are given in section 5. Benchmarks for computing Lie symmetries are given in 
section 6. We close the paper with some concluding remarks in section 7. 



2 Methods 

2. 1 Lie symmetries of differential equations 

Let {u\, . . . ,u n } = u be a set of functions (dependent variables) of the (in- 
dependent) variables {x\, . . . ,x m } = x. A system of p differential equations 
satisfied by the n functions Uj(x) can be written as 





(1) 



with 



Uji — 



d k Uj/dx h ■ ■ ■ dx ik , I = ii,..., z fc . 



(2) 



For m — 1 eq. (TO is a set of ODE's. 



A transformation of variables 




(3) 
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is a symmetry transformation of eq. ([T]) if 

F,(u' j ,x' t ,u' jI ) = 0, 



(4) 



where u'jj = d u'Jdx'^ ■ ■ -dx' ik , whenever eq. ([I]) holds, i. e. if eq. ([I]) is form 
invariant under Q, or equivalently, if eq. ([3]) maps a solution into another 
solution of eq. ([I]). Such transformations are called Lie symmetries (point 
symmetries). The set of Lie symmetries of a (system of) differential equa- 
tion^) is a Lie group, and therefore can be obtained from the knowledge of 
the infinitesimal transformations (in fact only the subgroup of transformations 
connected to the identity transformation) [10]. The infinitesimal symmetries 
can be written as: 



x\ = Xi + eOiiu, x), 

u'j = Uj + er]j(u,x), (5) 

where e is an infinitesimal parameter and 6*j and rjj are functions of the de- 
pendent and independent variables. The infinitesimal symmetry generator of 
transformation ([5]) is 

noma 

g = 2>^ + (6) 

j=l UU J i=l ^ 



The set of all infinitesimal symmetry generators form a Lie algebra with re- 
spect to the commutation operation. 

In order to determine the invariance condition of ([T]) under the infinitesimal 
transformation (151), we note that 



dx'i dx; 



+ e 



du i 86 1 



dxi 



i=i 



dxi dxi 



du 3 

dxi 



(i) 



(7) 



The transformation rules for higher order derivatives can be obtained similarly. 
In the general case we have: 



d k u' 3 _ d k Uj 
dx' h ■ ■ ■ dx' ik dx h ■ ■ ■ dx ik 



+ e Vj%-i k = Uj,ii~i h + enj%.4„, (8) 



with T)jJ 1 .,. iK functions of the independent and dependent variables and its 
derivatives [TU]. Supposing that the highest derivative in ([T]) is of order k we 
can express its invariance under an infinitesimal transformation by 

G (fe) F M = 0, (9) 
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where G^ fc ^ is the fc-th prolongation of the generator G in (jH]) and is obtained 
from eq. (jSJ) as: 

m Q n a k a 

1=1 j = l ""j |=lti,..,ii ua 0,lx-H 



Using the orthonomic form of eq. ([T]) to eliminate highest order derivatives 
from eq. (Q, and equating to zero the coefficients of the remaining derivatives, 
or more precisely the coefficients of linearly independent functions of the latter, 
we obtain the determining system for the symmetry transformations of eq. ([T]). 
The reduction to the orthonomic form of eq. (pQ) is performed using standard 
methods (see [42 43J and references therein). If the reduction is not possible 



SADE will issue an error message before aborting the calculations. 



The symmetry generator in eq. 
form: 



is equivalent to the following evolutionary 



g = E 



m dv ■ 



dxi 



_d_ 

8Ua 



i=i 



_d_ 

diij 
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Both forms as given in eqs. and (ITT]) describe the same transformation, in 
the sense that they map a given solution to the same transformed solution. 



2.1.1 Quasi- Polynomial Symmetries 

First order differential equations usually admit an infinite dimensional Lie 
symmetry group. To determine their Lie symmetries it is usually necessary to 
impose an ansatz on its symmetry generators. One possibility is to suppose 
that the coefficients 77, of the symmetry generator are polynomial functions of 
the dependent variables. A more general ansatz consists to consider the class of 
quasi-polynomial functions, as introduced in ref. [25], previously implemented 
by the authors in the package QPSI [24], and now included in SADE. 

A system of ODE's of the form 

m n 

±i = XiY,Aj Y[x k B ^ k ; i = l,...,n, (12) 

j=i k=i 



is called Quasi- Polynomial (QP) [32]. In eq. ( jl"2[) and Bj^ are real or com- 
plex constants and m is the number of different quasi-monomials in (I12p . We 
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define a new set of variables by the Quasi- Monomial Transformation (QMT): 

n „ 

yi=U x k Bik , (13) 



k=l 



where Bik — B^ for all % and k < n, Bik — for i < n and k > m, and 
Bik = $ik f° r n < i,k < m, in such way that the inverse transformation is 

also a QMT with the exponent matrix given by B (the case with B singular 
can also be handled as discussed in [32J). System ([TJ is cast by transformation 
(TT3]) into a quadratic system of equations, the Lotka-Volterra form: 

m 

yi = yi^2M ij y j} i = l,...,m, (14) 



where the matrix M is given by: 

M = BA. (15) 



Lotka-Volterra equations are extensively studied in the literature. Many results 
obtained for this special class of equations can then be recast into the more 
general QP form |25I27I50152155I53] . 

The central result obtained in Ref. [25] is that any quasi-polynomial symmetry 
generator G, such that 9^ — and rji a quasi-polynomial function of the 
dependent variables, can be decomposed as: 

G = ^G W , (16) 



with: 

G(*)=^ (i) TW, [F,G (i) ]=0, (17) 



where F is the flow associated to the quasi-polynomial system (fT2|) : 

m Q 

F = M ijyi y j7 r-, (18) 

i,j=i °yi 



y^ l) is a quasi-monomial: 

=(/...(^, (19) 
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with ^ real numbers (see eg. [221 below) and TW a polynomial semi-invariant 
vector field satisfying: 

[F,TW] = A W TW. (20) 

The eigenvalue X^ is a linear function of the form 

A W =EA?V (21) 



Both Aw and £W are solutions of the equation: 

EefM ife = -A«. (22) 



It is straightforward to show that if a symmetry generator (including the 
flow F) can be written as a linear combination of the remaining generators, 
then the coefficients of the expansion (as functions of Xi) are first-integrals of 
system fll2l) . 



2.2 Lie-Backlund symmetries 



Lie symmetries are diffeomorphisms on the space of dependent Uj and inde- 
pendent Xi variables. Lie-Backlund, or generalized, symmetries depend also on 
derivatives of Uj. We present here a brief account of how to compute the gener- 
ators of Lie-Backlund symmetries (for more details see [10J and [7J). The deter- 
mining equations for Lie-Backlund symmetries are more easily obtained using 
the evolutionary form (iTTj) . with rjj = r)j(iij,Xi,Uji) and 9{ = 9i(uj,Xi,Uji), 
with prolongation 

G W =EE^ ( 23 ) 



where I = . . . , i m , Dj = d k /dx^ ■ ■ ■ dxi k , Dq = 1, and the summation over / 
is taken for all values of indices such that \I\ — + . . .+ik < k. The invariance 
condition can be expressed as 

G^F^x,u, Ul )=0, (24) 



with k the maximum differentiation order of Uj in F^. In order to equate coef- 
ficients of independent derivatives of Uj in (Blj) , we distinguish dependent and 
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independent derivatives of Uj using the original system in the orthonomic form, 
and its differential consequences. In this way only independent derivatives re- 
main, and at this step each coefficient of derivatives that are not arguments of 
Qj is equated to zero, yielding the determining system for Lie-Backlund sym- 
metries. Its solution demands a greater computational effort than the solution 
of the analogous determining system for Lie symmetries. 

2.3 Reduction of PDE's and invariant solutions 

Symmetries of a differential system can be used to construct analytical solu- 
tions or a reduction into a system depending on a smaller number of indepen- 
dent variables. A symmetry generator as given in ([6]) can be transformed, by 
a change of dependent = Ti(u,x) and independent variables Sj = Sj(u,x), 
called canonical coordinates, into the form 

g, = A. (25 ) 

Solutions invariant under the symmetry generated by Gi do not depend on 
Si. In this way we obtain a reduction into a system with n — 1 independent 
variables. More generally p < m symmetry generators can de used to reduce 
to a new system with m — p independent variables, on the condition that 
a set of mutual canonical variables exists for the set of p generators. For 
p = m — 1 we obtain a system of ODE's. The latter, if solvable, then yields an 
analytical (particular) solution for the original system. In practice, considering 
p generators G^, % = 1, . . . ,p, we look for solutions Ui(x) satisfying 

G®uj(x) = 0, j = l,...,n, (26) 

where G^ is the evolutionary form (fTTT) . Equation (1261) is a linear system 
usually simpler to solve than the original system using the characteristics 
method. Replacing its solution into the original system yields a reduced system 
with n — p independent variables. 

2.4 Nonclassical symmetries 

Lie symmetries maps the set of all solutions of a differential system into itself. 
Invariant solutions then correspond, among all solutions, to those that are 
invariant under one or more symmetry transformations. More generally, non- 
classical symmetries transform a solution, still to be determined, into itself. 



9 



This amounts to require that both ([T]) and the invariance condition: 



Qi = G u,i(x) = rji-22 o-^j = 0, 



(27) 



j=i ux j 



are invariant under (jSJ). This is expressed by: 



G«F M = 0, 



(28) 



and 



G (fc) g, = o. 



(29) 



A more detailed account of nonclassical symmetries is given in Ref. [35] ■ It 
can be shown that eq. (129|) holds whenever eq. (128]) is satisfied. The result- 
ing system (128]) is non-linear in the unknowns rji and 8 { as it must be solved 
modulo eq. (129 p . and thus much harder to solve than the linear determining 
system for Lie symmetries. The set of all nonclassical symmetries include all 
Lie symmetries, and do not form a vector space (no associated Lie algebra). 
Another useful property is that if G is the generator of a nonclassical symme- 
try, then F(u, x) G also generates a nonclassical symmetry, for any arbitrary 
(sufficiently different iable) function F. As a consequence and without loss of 
generality, we consider the cases with 8i = 1 or 8\ = 0. In the later case there 
are still two possibilities: either 82 = 1 or 82 = 0, and so on. 

Computer algebra determination of the invariance condition for nonclassical 
symmetries can lead to infinite loops when replacing dependent derivatives 
from eq. (I2"9"j) into eq. (12S1) [36J . This is avoided in our approach by the following 

steps: 

(1) Chose an independent variable X}.. 

(2) Solve the invariance conditions (127)) for all derivatives, dui/dxk , % = 
l,...,n. 

(3) From the result of the previous step, eliminate all derivatives with respect 
to Xfc in ([1]). 

(4) Determine the invariance condition using the resulting differential system. 

(5) Replace in the determining system all derivatives of U{ with respect to Xk 
using step (J2J). 

The variable Xk is chosen such that MaxDer(xk) < MaxDer(xi) for i ^ k, 
with MaxDer(xi) the maximum derivative order of any dependent variable 
with respect to in ([T]). This usually results in a "simpler" determining system. 
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2.5 Potential symmetries 



For systems in conserved form, potential symmetries can be used to construct 
invariant solutions that are not obtainable neither from Lie nor nonclassical 
symmetries [6.37J. A partial differential equation is said to be in a conserved 
form if it can be written as: 

m q 

zl &r F i( u J> x i- u n) = °- ( 30 ) 
i=i ax * 

This implies that there exists m(m — 1) components (potentials) j (i < j) 
of an antisymmetric tensor such that 

m fl\Cr.. m /Mr.. 

f<= E(-i)^+ Ef-rJ. (si) 

i<j=l Udj i j<i=l 3 

The generalization to a system of PDE's is straightforward (each equation 
must be put in a conserved form). Eq. (13TT) is a system of M PDE's with 
1 + m(m — l)/2 dependent variables (v,j and the potentials). For M > 3 the 
system is under-determined, and some "gauge" conditions on the potentials 
must be given |37] . 

An infinitesimal symmetry of eq. (13TD 

u'i = Ui + erji(u, *S?,x), 
% 3 = % + e^(u,^,x), 

x'i =Xi + e9i(u,V,x), (32) 

is a potential symmetry of the original system (13"U|) if rji or $i depend on 
As a result the transformation for u and x« is non-local, since it depends on 
the potentials, which are solutions of eq. (13T|) . A potential symmetry can then 
be used to reduce eq. (130]) and, in some cases, to obtain invariant solutions. 



2.6 Symmetries and conservation laws 
2. 6. 1 Quasi-polynomial first-integrals 

For the special case of quasi-polynomial first-order systems, we first obtain 
the associated Lotka-Volterra form (Il4p . and define a semi-invariant (Darboux 
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polynomial) as a polynomial function f(y) such that 
m df 

/ = Z>«j£ = V. (33) 



where the eigenvalue A is a function of yi, . . . , y m . 

For the Lotka-Volterra form the following properties were proved in [25|33] : 

(i) A is a linear function, i.e., 

m 

A = £ X jVj . (34) 

(ii) Any polynomial semi-invariant / can be decomposed as: 

/ = £ / (P) , (35) 
p 

where is a homogeneous polynomial of degree p. Furthermore, each 

is also a semi-invariant with the same eigenvalue as /. Any Quasi-Polynomial 

first-integral can be decomposed as: 

J = y£jf»f<P), (36) 
p 

where y^ P f&> is a first-integral with a solution of eq. ( 122]) . 

If one of the Quasi- Monomials in fl5} is a constant, A may also admit a constant 
value with respect to the original variables Xk, and a first-integral can be 
obtained by multiplying the corresponding semi-invariant by exp(— At). It is 
easy to show that for f^ 1 ' and semi-invariants with respective eigenvalues 
A« and A (2) , /W/( 2 ) and fW/f® are also-semi-invariants with eigenvalues 
A* 1 ) + A^ and A^ — \( 2 \ respectively. The first integrals are then obtained 
by combinations of the form 

QP 1 {x)[QP 2 {x)f 1 exp{pt), (37) 



such that it has a vanishing eigenvalue, and therefore has a zero time deriva- 
tive. Analogously, it is straightforward to implement the computation of first- 
integrals of the form Pi(x)+log(x^) and P2(x, ln(x)), where Pi and P 2 are poly- 
nomials in their arguments. A similar result also holds for quasi-polynomial 
symmetries [24|. 
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2.6.2 Ncether theorem 



Many systems of interest are described by equations that can be deduced from 
a variational principle, with action S defined by 



S 



d m x. 



(38) 



where Uj (J = l,...,n) are the dependent variables, Xi (z = l,...,m) the 
independent variables and Ujj = duj/dxi. The lagrangian of the system is 
denoted by C and C is a bounded region of the m-dimensional space of inde- 
pendent variables. The differential system is obtained from the requirement 
that S is an extremum for any solution Uj(x). 

Nother theorem [38J states that every symmetry transformation of the action 
S of the form 



x ■ = %i + e9i{u,x), 

u'j = Uj + e rjj(u, x), (39) 
is related to a conservation law. The invariance of S under (1391) implies that 

VL = (40) 



with 

^ „ 9 x , d / dr)j x d6 k \ d ^ d9i , , „ . 

p s ? + ? + 5 (5 - ? ^ + ? (41) 



where /j are also unknowns to be determined from condition ( 140]) alongside 
with the 0j's and the r] M 's. The associated first-integral or conserved current is 
given by 

I i = £0 i + y £^L U - Y, «;,A ) - f u (42) 

j Ji* \ k / 



which satisfies the conservation law 
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2. 7 Reduction of order of an ODE 

Let us consider an ODE x^ = /(£, x, x', . . . , x^ k ~^) where x^ is the k-th 
derivative of x with respect to t, admitting a symmetry generator Gi. Using 
the canonical coordinates r (dependent variable) and s (independent vari- 
ables), we have Gi = d/dr, and consequently the original equation is cast in 
the form r^ = g(s,r', . . . , r^ -1 )), for some function g, which is an ODE of 
order k — 1 in u = r' . 

Now suppose the original equation admits another symmetry generator G2. 
It can also be rewritten using the same canonical variables. If the extended 

(k) 

generator G 2 is such that it does not act on s directly but only on its deriva- 
tives, then it can be used for a further reduction of order. This procedure can 
then be iterated for any number of generators, provided that at each step the 
generator used acts only in the remaining derivatives. Of course if m = n the 
system can be completely solved. This reduction is possible iff the Lie algebra 
spanned by the k generators is solvable |39| . 



2.8 Equations admitting a symmetry group 

In many situations the equations describing a given system are not known 
in closed form, but some of its symmetries are known. This is the case for 
instance if an underlying kinematical group (e. g. the Poincare or Lorenz 
group) is imposed by the physics of the system. This is frequently the case 
for transport equations for which no complete general theory exists |40j . One 
may hope that using the knowledge of symmetries may determine a class of 
equations for the problem at hand. In this way, Let us consider a set of p 
equations on the unknowns m, i — 1, . . . , n, of the form: 

f>X = 0, i = l,...,p, (44) 
i=i 



where Aj is a differential operator which can be non-linear. The system (jUJ) 
defines a class of equations if the operators depend on unknown functions. 
Now we impose that eq. (T4"4"|) is invariant under symmetry transformations 
generated by G«; i = 1, . . . , k, forming a system of differential equations for 
the unknown functions in A^, that can be solved in some cases. This is imple- 
mented in SADE in the routine equivalence. The main shortcoming here is 
that the system to be solved is non-linear. 
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3 Heuristics for solving the determining system 

No proved fully general, finite and terminating algorithm exists for the so- 
lution of linear or non-linear over-determined systems of partial differential 
equations, of the form obtained as determining systems for Lie, Lie-Backhand 
and nonclassical symmetries. Here we present the heuristics used in the present 
package. For linear systems, it was able to efficiently solve all test cases, span- 
ning a large number of equations ranging from the simple harmonic oscillator 
to equations describing coupled relativistic fields. 



3.1 Over- determined system of linear partial differential equations 

The basic idea is to solve simpler equations first, always trying to simplify 
further the system. Of course the meaning of "simpler" is quite subjective and 
our definition will become clear below. In a few cases it is necessary to append 
the determining system with integrability conditions for some, or all, of its 
equations. This is done using the MAPLE built-in routine rifsimp, when it is 
applicable, or otherwise using a slightly modified version of the Kolchin-Ritt 
algorithm with sorting [41 J. The reduction to the involutive form (see |4"2"f4"3"] 
for a proper definition) is usually very expensive in computational time, and 
should be done only if the system cannot be solved otherwise, and after some 
preliminary simplifications. Another strategy is to reduce only a subsystem 
with equations containing up to a prescribed number of terms. 

In what follows the number of terms in an equation is the number of its sum- 
mands. Parameters controlling the flow of the solution algorithm are specified 
in global variables that can be modified by the user, and with default values 
given below. These are the main steps of our algorithm: 

(1) Solve all algebraic equations with a maximum of 2 terms. 

(2) Solve all differential equations of the form 

d k f 

OXi l . . . %i k 

where / is any unknown in the determining system. 

(3) Solve any algebraic equation in the original unknowns 9i and r^. 

(4) Reduce to involutive form the subset of equations with at most A r 1 terms. 

(5) If any equation can be written as an expansion in linearly independent 
functions, then equate each coefficient to zero. Repeat until no more such 
decomposition is possible. 

(6) Repeat step [2j If any equations is solved, go to step dUl 
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(45) 



(7) Completely reduce to involutive form, and in case it succeeds, go to 
stepdS 

(8) Solve all equations with at most N 2 terms that can integrated as ODE's 
in one of the unknowns. If no equation can be solved, then repeat with 
A^2 + 3 terms, and so on up to the maximal value N%. If any equation is 
solved, then go to step [TUl 

(9) Solve one ODE with any number of terms. If it succeeds, go to step [TJ 

(10) Repeat step [I] 

(11) Look for all equations that are expansions in linearly independent func- 
tions, and equate to zero each coefficient of the expansion. 

(12) Repeat step [I] 

(13) Repeat step[2j 

(14) Repeat stepEJ 

(15) Solve all algebraic equations (with any number of terms). If any equation 
is solved, then go to step [JUJ 

(16) Reduce to the involutive form. If succeeds, go to to step [TOl 

(17) Look for one equation of the form fx(xi,x 2 ) = fi{xii%?i) and replace it 
by fi(x!,x 2 ) = h{x\) and f 2 (x 1 ,x 3 ) = f 3 (x 1 ), where f 1 and f 2 are two 
unknowns in the system and f 3 a new unknown. If succeeds go to step [TIJJ 

(18) Repeat step El 

The whole algorithm is repeated until the system is completely solved or no 
additional simplification occurs. The following global variables are used: 

Ni = SADE[partial_reduction] 

N 2 = SADE[_ne] 

N 3 = SADE[_nmaxeq] 

with default values N% = N 2 = 5 and N 3 = 8. 



3.2 Non-linear systems 

Solving non-linear overdetermined systems of PDE's, as those resulting from 
the determination of nonclassical symmetries, is a very difficult task. Its imple- 
mentation in SADE is still under development, but can nevertheless be used 
in some interesting cases. There are other more efficient algorithms, such as 
the one used in the REDUCE package CRACK [13J. Our algorithm can be 
roughly sketched as: 

(1) Solve all linear equations. 

(2) Reduce the resulting system to the involutive form. 

(3) Solve all linear equations. 
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(4) Decomposes equations which are expansions in linearly independent func- 
tions. 

(5) Solve a single nonlinear ODE. If any equation is solved, repeat step H] and 
go to step [TJ 

(6) Solve all purely algebraic equations, considering multiple solutions. If any 
equation is solved, then go to step [TJ 

There are also options for reducing the determining system to involutive form 
before trying to solve it, and to use the MAPLE builtin routine for solving 
PDE's. 



4 Package Commands 

Here we briefly describe each command available in SADE. The examples given 
in section [5] should be self explanatory and complementary to this section. The 
following abbreviations are used for the input arguments: 

eqs: a single or a set of differential equations, 
unks: the unknowns in eqs. 

gen: a symmetry generator, written in SADE notation (see section [5]) . 

depvars: list of dependent variables. 

indepvars: list of independent variables. 

vars: list of dependent and independent variables. 

der_order: list with the derivative order of each dependent variable in the 

independent variables (see section [578]) . 
drvs: set with derivatives of unknowns in eqs. 
subs _ rule: a substitution rule, 
funcs: a set with new undetermined functions, 
name: a maple variable name, 
opt: optional arguments. 

determining: optional argument to return only the determining system; 
involutive: optional argument to reduce the determining system to involu- 
tive form, 
params: a set on free parameters. 

Package commands and corresponding inputs are given in the following listing: 

liesymmetries(eqs, unks, opt): Computes Lie symmetry generators. Op- 
tional arguments: determining, involutive, case= n, n integer, only the 
case with 9 n = 1, 9i = (i < n) is considered, builtin - solves the deter- 
mining system using the MAPLE builtin command pdsolve. 
def ault_parameters - the determining system is solved using default pa- 
rameters reducing CPU time, although the system may not be completely 
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solved. 



ncsymmetries(eqs,unks,opts): Computes nonclassical symmetry genera- 
tors of DE's. Optional arguments: determining and involutive. 

LBsymmetries(eqs,unks,opts): Obtains Lie-Backlund symmetry genera- 
tors. Optional arguments: determining, involutive and 
parameter = paramset - computes the generators with conditions on the 
free parameters in paramset. 

lindsolve(eqs,unks): Solves a linear overdetermined system of PDE's. 

nonlindsolve(eqs,unks): Solves a non-linear overdetermined system of par- 
tial differential equations. 

casimir_ invariant ({geni, gen 2 , . . .},depvars,indepvars,der_order): 

Computes the Casimir invariants of a set of generators. 

ansatz(subs_rule,funcs): Applies a set of ansatze to the determining equa- 
tions. This routine can be used either once the determining system is ob- 
tained or if SADE could not completely solve the determining system. The 
elements of subs_rule must be given in the form function = expression 
with new undetermined functions in expression specified in funcs. 

noether(lagrangian,funcs,gen): Computes the Nother conserved currents 
or first-integrals from a lagrangian function. 

equivalence(eqs,{geni, gen2, . . .}, funcs): Obtains the most generic form 
of a class of equations admitting a symmetry algebra. 

comm(gen 1 , gen 2 ,vars): Commutator of two linear operators (generators). 

com_table({geni, gen2, . . .},vars,name): Commutation table of a set of 
infinitesimal generators, name specifies a name to represent each generator. 

AdjointRep({gen!, gen 2 . . .},vars, name, parameter): Computes the ta- 
ble with the action of adjoint maps on each generator of a Lie Algebra, name 
is used to represent each generator in the table and parameter specifies the 
adjoint Lie group parameter 

StructConst({gen!, gen 2 . . .},vars): Computes an array with the struc- 
ture constants of a Lie algebra. 

linear_rep(operator,{geni, gen 2 . . .},vars,name): Determines the most 
general linear operator representing a class of differential equations defined 
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by operator admitting a symmetry algebra. 

PDEreduction(eqs,unks,{gen 1 , gen 2 •••}): Obtains the reduced form of 
a PDE or a PDE system from a set of symmetry generators. For a sys- 
tem with M independent variables K symmetry generators can be used to 
reduce to a new system depends with M — K (transformed) independent 
variables. 

invar iant _sol ( eqs, unks,{ geri!, gen 2 .. .}): Obtains invariant solutions of 
a PDE or a system of PDE's using symmetry generators. 

issolvable({geni, gen 2 . . .},vars): Tests if a Lie algebra is solvable. 

canonical_basis({gen!, gen 2 . . .},vars): Computes the canonical basis of 
a Lie algebra. 

derived _ subalg({ geni, gen 2 .. .},vars): Computes the generators of the 
derived subalgebra of a Lie algebra. 

odesolver(eqs,{gen!, gen 2 . . .},unks,opt): Solves an ODE by successive 
reductions using a solvable Lie algebra. Optional argument: transformation 
- returns only the transformation of variables solving the system. 

reduce_ode_sist(eqs,{gen!, gen 2 . . .},depvars 1 ,indepvar 1 ,depvars 2 , 
indepvar 2 ): Reduces by one the dimension of a system of first order ODE's 
using a symmetry generator. Note that new dependent and independent 
variables must be given and are represented by the index 2. Index 1 denotes 
original variables. 

ode_reduce_orderl(eqs,gen,depvar 1 ,mdepvar 1 ,depvar2,mdepvar 2 ): 

Reduces by one the order a a single ODE using a symmetry generator. 

ode_invsolution(eqs,unks,gen): Obtains invariant solutions for a single 
ODE. 

conserved(eqs,unks,params,n,opt): Obtains the QP-invariants of a QP 
first order system by computing a Darboux polynomial up to degree n. Op- 
tional arguments: Groebner - a Grobner basis computation is used to solve 
the polynomial system of determining equations, positive - the results are 
simplified to the positive orthant. surfaces - returns the defining equations 
for invariant hyper-surfaces. 

QPsymmetries(eqs,unks,params,n): Determines QP symmetry genera- 
tors with n the degree of the polynomial quasi-symmetry |24] . 
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verif_if_inv(eqs,unks,params): Determines parameter values such that 
non-trivial QP first-integrals (i. e. with non-integer exponents) may exist. 

constlog(eqs,unks,params,n): Computes first-integrals of the form P 1 (x)+ 
log(a^) and P 2 (x,ln(x)), with Pi and P 2 polynomials of degree n and x 
stands for all dependent variables. 



5 Illustrative examples 

We present some illustrative cases of a basic use of SADE. A whole suite of 
examples, including more complex problems, is given with the distribution 
files. 



5.1 Lie symmetries 

Let us consider first as a simple example the heat diffusion equation: du/dt = 
d 2 u/dx 2 . The input to determine its Lie symmetries and the corresponding 
output are: 

> g := liesymmetries({dif f (u(x, t), t) = dif f (u(x, t), x, x)}, [u(x, t)]); 

g := [{D t , D x , uD u , _F1 (t, x) D u , 2 tD t + xD x , uxD u - 2 tD x , 
t 2 D t + 1/4 (-2ut - mx 2 ) D u + txD x } , 

-|- w <'-4 

The first element of the list in the output is the set of infinitesimal symmetry 
generators. The second element is a list of constraints on these generators, 
which in this case is the original equation. The element D[a] — > D a in the out- 
put stands for the directional derivative 8/ da, where a is any of the dependent 
or independent variables. The commutator of two generators is computed as: 

> comm(D[x], — 1/2 * u * x * D[u] + t * D[x], [u, x, t]); 




The complete commutation relations for the finite dimensional algebra can be 
computed as follows: 

> gen := convert(g[l] minus {_Fl(t, x) * D[u]}, list); 
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> com_table(gen, [u, x, t],G); 
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5.# Lie-Backlund symmetries 

To illustrate how dependencies on derivatives of dependent variables are han- 
dled when computing Lie-Backlund symmetries, we consider the following two- 
dimensional PDE system [9]: 

d 2 u 2 _ ldu 2 d 2 ui _ dux u\ 
dx\ 2 dx 2 ' dx\ dx 2 2 ' 

Requiring, for instance, that the evolutionary form of the symmetry generators 
depends on dui/dxi, d 2 U\jdx\, du 2 jdx\ and d 3 u 2 /dxl, implies that the in- 
dependent derivatives are dui/dxi, dui/dx 2 , du 2 jdx\ and d 2 u 2 jdx\dx 2 . Here 
are the corresponding input and output: 



> LBsymmetries(eqs, [ul(xl,x2),u2(xl,x2)], 

{diff{ul{xl, x2), xl), diff(ul(xl, x2), xl$2), diff(u2(xl, x2), xl), 
diff(u2(xl,x2),xl$3)}); 



1/2 
+1/2 



(xl,x2) + 2_Fl (xl,x2)u2)D ul 
u2 
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[2^_Fl( X l, X 2)-^- 2 _Fl( X l, X 2) 

1 dxl ' dx2 ' dxl ' <9a;l <9z,2 

The first element is a set with the Lie-Backlund symmetry generators, the 
second element a set with constraints on the generators (an empty set if none), 
and the third element the set of independent derivatives. 



5.3 Solving an ODE 



Let us consider the equation [SJ: 

fdu\ 5 d 3 u ^ fdu\ 4 (d 2 u\ 2 fd 2 u\ 3 ^ 
\dx J dx 3 \dx I \ dx 2 J \ dx 2 J 



with u = u(x). First compute its Lie symmetries: 

> eq := diff(u(x),x)"5 * diff(u(x), x, x, x) — diff(u(x),x,x)~3 

-3* diff(u(x),x)"4* diff(u(x),x, x)"2 : 
liesymmetries(eq, [u(x)]); 

[{D u ,D x ,uD x ,uD u + 3/2xD x },{}] 
then chose a solvable subalgebra: 

> ls:= {D u ,D x ,uD x } : 
issolvable(ls, [u,x]); 



true 

Here are the inputs for computing the solutions and displaying the first one 
(there are three branches): 

> sol := odesolver(eq, {D[u], D[x], u * D[x]}, [u(x)]) : 
>sol[l]; 

u{x) = 1/2 _C1 - 1/2 (1/2 (-12x + 6 _C2 _C1 - 12 _C3 - _C2 3 

+2 (36 x 2 - 36 x_ C2 _ CI + 72 x _C3 + 6 x_C2 3 

+9 C2 2 C1 2 -3Q C2 CI C3-3 C2 A (71+36 C3 2 
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+6 _ C3 _ C2 3 ) 1/2 ) 1/3 + 1/2 _ C2 2 /((-12 x + 6 _ C2 _ CI - 12 _ C3 
-_C2 3 + 2(36x 2 -36x_C2 _C1 + 72 x_C3 + 6 x _C2 3 
+9 _C2 2 _C1 2 - 36 _C2 _C1 _C3 -3 _C2 4 _C1 + 36_C? 2 
+6 _ C3 _ C2 3 ) 1/2 ) 1/3 ) 1 /2 ) 2 

5.-^ Invariant solutions 



The Burgers equation in one dimension can be rewritten as a set of two one- 
dimensional equations by defining v = du/dx |44j: 

9u du dv 

__ B = ; =0. (48) 



Since the number of independent variables is two, only a single symmetry 
generator is necessary to obtain a group invariant solution. In this form, one 
of the symmetry generators admitted by the Burgers equation is: 

G — (ut- x) D u + {2vt- 1) D v - t 2 D t - txD x . (49) 



The associated invariant solution is obtained using the following input in 
SADE: 



> invariant_sol(eq, [u(x, t), v(x, t)], {(u * t — x) * D[u] 
+ (2 * v * t - 1) * D[v] - t 2 * D[t] - t * x * D[x]}); 



{{u{x,t) = -{-_C1 x + tanh(l/2 ~ ))_C1 

Cj 1 t 

T _l no f 

v(x,t) = 1/2 (2t(cosh(l/2 X + l )f_Cl 2 - l)r 2 
(cosh(l/2^±^))- 2 _CT- 2 }} 



5.5 PDE reduction 



Let us consider a massless nonlinear Klein-Gordon equation with a A0 4 /4 
self-interaction potential: 
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> p := phi(x,y,z,t) : 

eq := diff(p, x, x) + diff (p, y, y) + diff (p, z, z) — diff(p, t, t) + lambda * p 3 : 

This equation can be reduced to a PDE with two independent variables using 
the symmetry generators: 

G x =yD t + tD y , G 2 = zD x -xD z , (50) 

a Lorenz boost in the y direction and a spatial rotation around the y axis, 
respectively. The reduced equation is obtained using the SADE command: 

> PDEreduction(eq, [phi(x, y, z, t)], yD t + tD y , zD x — xD z ); 



{[{(f>(x,y,z,t) = _Fl(x 2 + z 2 ,-y 2 + t 2 )}, 

{ 4 (^j2- F1 (^^2mi+4^j_F1^1^2)-4^_F1^1^2) 

"4 (^2 -Fl (tl , + A (_F1 £1 , £2)) 3 }, 

{£1 = x 2 + z 2 ,£2 = -y 2 + t 2 }]} 

The output is a set with the different possible reductions (one in the present 
case). Each element is a list: the first element defines the relation between 
the original unknown(s) and the reduced dependent variable(s). The second 
element is the set of reduced equations, and the last element is a set with the 
similarity variables £j used to reduce the original equation(s). 



5.6 Nother theorem 



As an example, we consider the relativistic massless scalar field in 1+1 dimen- 
sions, with lagrangian density: 

c = w-w = - (51) 



In this case a conserved current I is a two-vector function of 0, x and t such 
that dli/dx + dlijdt = 0. They can be obtained in the following way: 

> noether(diff (phi(x,t),x)~2/2 - dif f (phi(x, t), t)~2/2, [phi(x, t)]); 
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'd_\d_ 



n 9 n d 
—2 —ip, 2 —ip 
at ax 



d \ d H ln ( d ^ 



dt dx 



dx 



> [t,x] 



Here the output is a sequence. The first element is the set of conserved currents, 
and the second element the independent variables defining the component 
ordering used for the components of I. 



5.7 Equations admitting a symmetry algebra 



Let us consider the following family of PDE's: 

du + + d\du 2 ) = ^ 

dt dx dx 2 



where / = f(x,t) and al = d(x,t) are functions to be determined. Requiring 
that eq. ([52} admits G x = d/dt, G 2 = d/dx and G 3 = (x/2)d/dx + td/dt as 
symmetry generators restricts the postulated generic form (152]) : 

> gen := {D[t], D[x], (1/2) * x * D[x] + t* D[t]} : 

> eq := difF(u(x, t), t) + diff(f(x) * u(x,t),x) + diff(d(x) * u(x,t)"2,x,x) : 

> equivalence(eq, [u(x, t)], gen, {f (x, t), d(x, t)}); 



{[^u(x, t) + 2_Cl (^u(x, t)f + 2_C1 u(x, t)-^u(x, t), {}], 
[^u(x, t) + (^-f(x))u(x, t) + f(x)-^u(x, t), {}}} 

The output is a set of lists. In each list, the first element is a restricted form 
for the class of equations, the second element is the set of remaining equations 
still to be solved, none in this case. 



5.8 Casimir invariants 



The order of the derivatives in a Casimir invariant is specified by a list such 
that each element is a list with the maximum derivative order of the corre- 
sponding dependet variable in each independent variable. Note that the order 
of the derivatives in the invariants are one order higher than specified in the 
input (this is due to the specific algorithm used in the computation) : 
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Table 1 

Benchmarks for Lie Symmetry determination. 



System 


N eq 


SADE 


DESOLV 


PDEtools 


Heat equation |10) 


9 


0.4 


0.13 


0.4 


Klein- Gordon |45j 


32 


1.1 


0.69 


0.9 


Magneto-Hydro-Dynamics |45| 


39 


14.4 


12.9 


13.4 


Navier-Stokes [6] 


125 


3.9 


4.8 


2.9 


Dirac [9] 


352 


14.7 


20.1 




Gross-Neveu |46| 


352 


15.4 


211 




Maxwell-Dirac [9] 


2,621 


27.8 


1,121. 




Yang-Mills SU(2) 07] 


12,361 


82 


188 




Yang-Mills SU(3) g8] 


175,042 


13,847 







> casimir _invariant(x * D[t] + t * D [x] , v * D [v] + u* D[u] — x * D[x], 
HM, [[2,2], [2,2]]); 



(£>u(x,t)) (g u(x,t)) - (gf^ u(x,t) 
u (x, t) 4 

d d 9 3 d 3 

-u(x,t),-u(x,t), — u(x,t),— u(x,t), 

2xt (gf^ u{x,t)) +x 2 (£, u {x,t))+t 2 ^u{x,t))' 

u (x, t) 



6 Lie Symmetries Benchmarks 



It is beyond the scope of the present work to do an exhaustive comparison 
of SADE to other similar packages (see References [22,23] for an assessment 
for previous packages). We present only a comparison for the computation of 
Lie symmetries (the core of SADE), and leave a more thorough comparison 
to be presented in a future work. Since the release 11 of MAPLE it became 
possible to compute Lie symmetries using its native package PDEtools. The 
package GEM relies on the solution of the determining system using MAPLE 
built-in routines, leading to large computational effort in CPU time and mem- 
ory when solving the determining system, becoming intractable for systems 
such as Maxwell-Dirac or SU(2) and SU(3) Yang-Mills field equations. In this 
way, and considering only MAPLE packages, tabled] shows CPU times for the 
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determination of Lie symmetries (solving the determining system and obtain 
the symmetry generators in explicit form) of some representative systems us- 
ing SADE, DESOLV and PDEtools. All computations were performed on a 
i5 2.40 GHz computer, and using the corresponding automated routines for 
computing Lie symmetries in each package: symmetry and genvec in DESOLV, 
liesymmetries in SADE and Infinitesimals in PDEtools. The absence of a value 
in the table means that either the computation was not completed after a very 
long time, or that the package was unable to obtain all symmetries of the equa- 
tion. From table [T] we see that PDEtools handles only simpler systems (the 
three first lines in the table), with DESOLV being the fastest and PDEtools 
the slowest. Nevertheless, for the remaining cases in table [H with a number 
of equations in the determining system ranging from a few hundreds to hun- 
dreds of thousands, SADE performs better than DESOLV. For the Yang-Mills 
SU(3) field equations, SADE took a little less than 6 hours while DESOLV 
produced no output after more than 60 hours. On its turn, PDEtools missed 
some symmetries of the Navier-Stokes and all of the Dirac equations, and for 
the remaining equations was unable to terminate the computations after a 
considerable long time. We used MAPLE 14 for timings, except for DESOLV 
which performs better in MALE 13 (for the latest version available to the 
authors) . 



7 Concluding Remarks and Perspectives 



The present package implements symmetry methods for differential equations 
in MAPLE, including Lie, nonclassical, Lie-Backlund and potential symme- 
tries, the Quasi-Polynomial formalism and the computation of invariant solu- 
tions and reduction of ODE's and PDE's. SADE is well suited to handle more 
complicated systems, such as Maxwell-Dirac and Yang-Mills SU(2) and SU(3) 
equations, for which optimization is crucial. It also obtains all symmetries for 
the "difficult" Vaidya and Jacob- Jones systems |22| . We performed no direct 
performance comparisons besides those for Lie symmetries determination. The 
implementation of Markus algebras [34 25J as a tool for the computation of 
quasi-polynomial first-integrals, symmetries and invariant hyper-surfaces, and 
a better algorithm for solving overdetermined systems of non-linear PDE's are 
currently under implementation. We also hope to implement a more complete 
set of routines for the computation of conservation laws for PDE's. 
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